On phases in weakly interacting finite Bose systems 
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We study precursors of thermal phase transitions in finite 
systems of interacting Bose gases. For weakly repulsive in- 
teractions there is a phase transition to the one-vortex state. 
The distribution of zeros of the partition function indicates 
that this transition is first order, and the precursors of the 
phase transition are already displayed in systems of a few 
dozen bosons. Systems of this size do not exhibit new phases 
as more vortices are added to the system. 

PACS numbers: 03.75.Fi, 05.30.Jp, 67.40. Db 



I. INTRODUCTION 

The study of vortices in dilute atomic Bose-Einstein 
condensates has received a lot of attention in recent years. 
The condensate wave function of a vortex state exhibits a 
quantized circulation of its velocity field. This state may 
experimentally be formed, e.g., by "phase imprinting" 
techniques pj or by directly transferring angular momen- 
tum to the condensating system |2| . Current experiments 
are performed in the Thomas-Fermi regime of short co- 
herence length. For a review of recent results see, e.g., 
H . The regime of long coherence length is of interest as 
well. Indeed, harmonically trapped Bose systems with 
perturbatively weak repulsive interactions display a rich 
structure |Q]. At low ratios of angular momentum L to 
particle number N, i.e. L/N <C 1, the ground states are 
dominated by quadrupole and octupole excitations HU . 
The ground state structure changes smoothly until the 
one-vortex state L/N = 1 is approached. In particu- 
lar, the ground state wave function is known analytically 
and the ground state energy depends linearly on L for 
L/N < 1 ||. The one- vortex state is a Bose-Einstein 
condensed state where a macroscopic number of parti- 
cles carries one quantum of angular momentum each [H . 
The structure of ground states changes considerably as 
further angular momentum is put into the system. Re- 
sults from mean-field calculations exhibit changing sym- 
metries of the developing vortex array which correspond 
to changes in the occupation numbers of single-particle 
states M . This goes along with an increasing complexity 
of higher vortex states. The two-vortex state L/N = 2, 
for instance, already exhibits a complicated structure and 
is dominated by excitations that carry 0, 2, and 4 units 
of angular momentum H . Naturally, the question arises 
whether the observed changes in the condensate wave 
function structure and the formation of quantized vor- 
tices are associated with thermal phase transitions. 

For the transition to the one-vortex state, the answer 



is affirmative. Wilkin et al. M showed that the one- 
particle reduced density matrix of the one-vortex state 
has one eigenvalue of order N and thus meets a criterion 
for Bose-Einstein condensation in finite systems; see e.g. 
[|9|. The order of this phase transition, however, is not 
known until now. It is the main purpose of the present 
paper to classify this phase transition from its precursors 
in finite systems. 

We recall that small systems do not exhibit phase tran- 
sitions. Nevertheless, finite systems may display pre- 
cursors of phase transitions. Recently, Borrmann et al. 
proposed that the distributions of zeros of the canonical 
partition function "reveal the thermodynamic secrets of 
small systems in a distinct manner" [|10|. We will ap- 
ply this method to the case of weakly interacting Bose 
systems under rotation. 



II. SYSTEM AND METHOD 

We consider a system of N bosons confined in a three- 
dimensional harmonic trap at total angular momentum L 
and restrict ourselves to the sector of maximal magnetic 
quantum number, i.e. the particles are in the ground 
state with respect to excitations along the axis of rota- 
tion. The non-interacting system is highly degenerate. In 
what follows, we assume that the repulsive interaction be- 
tween the bosons is perturbatively weak and simply lifts 
this degeneracy. This yields sets of now quasi-degenerate 
states that are separated by multiples of the oscillator 
spacing Hu>. Under these conditions, the level spacing be- 
tween quasi-degenerate states is much smaller than the 
oscillator spacing. We are interested in the thermal prop- 
erties of the system while keeping the angular momentum 
L fixed. For temperatures that are smaller than the os- 
cillator spacing, we may restrict ourselves to the lowest- 
lying set of quasi-degenerate states with approximate en- 
ergies E k, LfiLO. (We set the ground state energy of the 
non-rotating system to zero.) Current experiments do 
not work within this low-temperature regime. We recall 
that the onset of Bose-Einstein condensation is already 
observed for temperatures kT ~ N x ^huj JLlJ. Below we 
find that quantitative results concerning thermal phase 
transitions are already determined by a few hundred low- 
lying levels. It thus seems that one can lift the require- 
ment of perturbatively small temperatures without fac- 
ing the inclusion of highly excited states for the problem 
under consideration. In the low-temperature regime, the 
Hamiltonian of the A^-boson system with contact inter- 
action reads M 
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Here, uo denotes the strength of the contact interaction. 
The operators dj and a,j create and annihilate one bo- 
son in the single-particle state with angular momentum 
j with j = 0, 1, 2 . . ., respectively. The operators 



N = 



£' 



and 



E 



j«j 



(2) 



(3) 



where n^ de- 



count the number of particles and the quanta of angular 
momentum and have quantum numbers N and L, respec 
tively. We have used the number operators hj 
Basis states are denoted as \no, n\, n%, . . .} 
notes the number of particles with angular momentum j 
Hilbert space is partition space, and the number of basis 
states is given by the number of partitions of L into at 
most N integers. The dimension of Hilbert space grows 
exponentially with increasing L while being only mildly 
dependent on N for L > N. 

Let us now turn to a description of the method pro- 
posed in Ref. [[l0|]l2[ |. It is based on an analysis that 
Grossmann and Rosenhauer made for macroscopic sys- 
tems about three decades ago |13|| . The canonical par- 
tition function Z(B) is evaluated at complex arguments 
B = (3 + it. Since Z(f3) is real, it suffices to consider 
the partition function Z(B) for arguments in the upper 
complex plane. Two different phases of macroscopic sys- 
tems are separated by a line of zeros which intersects the 
real axis at the critical temperature. Further information 
about the order of the phase transition is encoded in the 
slope of the line at the intersection point and the density 
of zeros close to the intersection point. This is physically 
plausible. We recall that thermodynamic quantities are 
given by logarithmic derivatives of the partition function 
and thus diverge at its zeros. 

In finite systems, the line of zeros reduces to more or 
less closely spaced zeros that line up on a curve. Fol- 
lowing Ref. [Hj, one then studies the behavior of the 
zeros with smallest imaginary part and of the underlying 
curve while increasing the number of particles. This al- 
lows one to predict the critical temperature and order of 
the phase transition in the infinite system. Furthermore, 
the shape of the curve or the presence of several such 
curves allows one to identify different "phases" already in 
finite systems. These techniques have been used to study 
precursors of phase transitions in Bose-Einstein conden- 
sates of ideal gases and atomic cluster jltjjl^]. Further 
applications include the classification of a pairing phase 
transition in finite Fermi systems J14| . 

We compute the canonical partition function from the 
eigenvalues of the Hamiltonian (||). This requires the 



complete diagonalization of H and is practicable only for 
moderately large values of L. Zeros of the partition func- 
tion are determined conveniently by locating the poles of 
the specific heat, which is given by 
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For larger system sizes (i.e. larger values of angular mo- 
mentum L), we restrict the computation to the lowest 
hundreds of levels and construct an approximate parti- 
tion function only. Nevertheless, we find sufficiently well- 
converged zeros of the partition function close to the real 
axis of complex temperature B = (3 + ir. This allows us 
to present quantitative data concerning the order of the 
phase transition in the macroscopic system. 



III. RESULTS 

We investigate the transition to the one-vortex state 
first. To this purpose, we fix the number of parti- 
cles to N = 30 and fully diagonalize the Hamiltonian 
for several values of angular momentum in the range 
0.8 < L/N < 1.07 around the one- vortex state L/N = 1. 
The partition function is computed from the obtained en- 
ergy levels. Figure [j] shows the specific heat as a function 
of complex B = f3 + ir. We restricted the plot to include 
only those zeros with smallest positive imaginary part. 
Figure [0 demonstrates that the zeros approach the real 
axis with increasing L. The closest encounter is found 
for the one- vortex state L = TV = 30. This is a precursor 
of the condensation into the one- vortex state in the infi- 
nite system and supports earlier results ||[7|||- We are 
particularly interested in the order of this phase transi- 
tion. To this purpose, we consider the one-vortex state 
L = N and compute eigenvalues of the Hamiltonian (jl]) 
for increasing values of particle number N. A complete 
diagonalization is prohibitively expensive for N exceed- 
ing values of about 35. Instead, we restrict ourselves to 
the computation of the lowest-lying eigenvalues. These 
are used for an approximate construction of the parti- 
tion function. We found numerically that its zeros with 
smallest positive imaginary parts are already sufficiently 
well converged when only a few hundred eigenvalues are 
included in the computation. We considered systems 
up to L = N = 55, corresponding to a dimension of 
Hilbert space of the order 4.5 x 10 5 . The relevant eigen- 
values are computed numerically using the ARPACK and 
PARPACK routines |lj|. 

Figure |2| shows the distribution of poles in the com- 
plex plane of the specific heat at the one- vortex state for 
different system sizes. These plots are generated from 
the lowest-lying 300 eigenvalues. (Increasing the num- 
ber of eigenvalues from 300 to 380 yields less than 1% 
change in the numerical results. Thus, the data is suffi- 
ciently well converged for our purposes.) It is clearly seen 
that the zeros line up and approach the real axis with in- 
creasing system size. The order of the phase transition 



is determined as follows [|T(J. The distribution of zeros 
close to the real axis is approximately described by three 
parameters. Two of these parameters reflect the order of 
the phase transition, while the third indicates the size of 
the system. Let us assume that the zeros lie on a line. 
We label the zeros according to their closeness to the real 
axis. Thus t\ reflects the discreteness of the system. The 
density of zeros for a given r k is given by 
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with k = 2,3,4, •■• A simple power law describes the 
density of zeros for small r, namely </>(r) ~t°. If we use 
only the first three zeros, then a is given by 

Q = ln0(T 3 )-ln(ft(T 2 ) 
In r 3 - In r 2 

The final parameter that describes the distribution of 
zeros is given by 7 = tan^ ~ (fa — /3i)/(t2 — T\). 

In the thermodynamic limit, 77 — > 0, in which case 
the parameters a and 7 coincide with the infinite system 
limits discussed by Grossmann and Rosenhauer J13| . For 
infinite systems, a = and 7 = indicates a first-order 
phase transition, while < a < 1 and 7 = or 7 ^ in- 
dicates a second-order transition. For systems approach- 
ing infinite particle number, a cannot be smaller than 
zero since this causes a divergence of the internal energy. 
In small systems, with finite 17, a < is possible and is 
also indicative of a first-order transition. We show our 
results for a, 7, 77 in Table | for the N — 40, 50, 52, 55 sys- 
tems. 7~i decreases with increasing system size as a power 
law. The fit is given by 17 = 1.15iV~ 153 . From the table 
we note that 7 is nearly zero, and a is a small negative 
number for each system we studied here. The critical 
temperature kT c is approximately given by l//3j.. Our 
data suggests that kT c ~ ugiV 11 - 1 - 4 . The iV-dependence 
differs considerably from what is found for the ideal Bose 
gas in three-dimensional traps [nTl. 

We thus find that the phase transition to the one- 
vortex state is first order. This is the main result of this 
work and combines well with previous results found for 
the energy and wave function structure. We recall that 
the ground state energy of the iV-boson system exhibits 
a kink at L = N, and that the wave function structure 
changes strongly when increasing L beyond N [JTjD . 

Are there further indications of phase transitions as 
one approaches the two-vortex state? Mean-field re- 
sults by Kavoulakis et al. p] show that the structure of 
ground states changes considerably at certain ratios of 
L/N. At L/N w 1.75 the mean-field density develops 
a two-fold symmetry, and the expectation rij for occu- 
pation of single-particle orbitals with odd j vanishes at 
L/N rs 1.75. A second change in ground state struc- 
ture occurs at L/N rs 2.03 where the two- fold symme- 
try of the mean-field density changes into a three-fold 
symmetry. This is accompanied by a macroscopic oc- 
cupation of single-particle orbitals with angular momen- 
tum j = 0, 3, 6, 9. We numerically computed the ground 



states of Hamiltonian (J]]) for systems with moderate par- 
ticle numbers and confirm these results. We also com- 
puted the partition function from the fully diagonalized 
Hamiltonian for L = 30 and N ranging from 10 to 20. 
However, the distribution of zeros of the partition func- 
tion does not indicate a thermal phase transition. In par- 
ticular, we do not find zeros that start to line up and ap- 
proach the real axis at the corresponding values of L/N . 
This is most likely due to the small system size. 

Let us finally turn to Bose-Einstein condensates with 
attractive interaction. In this case, the ground state is 
the (Bose-Einstein condensed) state of the non-rotating 
system which executes pure center-of-mass motion, i.e. 
all angular momentum is carried by the motion of the 
center of mass ]8|||. A direct calculation of the one- 
particle reduced density matrix indicates the presence 
of a condensate that is fragmented over several states 
H. Only when the center-of-mass motion is separated 
does one find the condensation into one state p6[ . The 
distribution of zeros of the partition function does not 
indicate a phase transition. This remains the case also 
when considering systems of size L, N rj 50. 

In summary, we have investigated phases in weakly 
interacting Bose-Einstein condensates. The phase tran- 
sition to the one-vortex state is first order, and its pre- 
cursors are clearly seen in the distribution of zeros for 
systems comprising only 30 to 55 particles. The question 
about further phase transitions to states with several vor- 
tices remains open. 
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N ti 7 a 

40 0.0715 -0.03 -0.1109 

50 0.055 0.0 -0.145 

52 0.053 0.003 -0.147 

55 0.050 0.003 -0.148 

TABLE I. Calculated parameters n, 7, and a for the various boson systems at L/iV = 1. 
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FIG. 1. Contour plots of the specific heat in the complex temperature plane for the N = 30 system at a) L = 20, b) L — 25, 
and c) L — 30 units of angular momentum. The spots indicate the locations of the zeros of the canonical partition function. 
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FIG. 2. Contour plots of the specific heat in the complex temperature plane for the L/N = 1 systems with increasing 
numbers of particles: a) N — 40, b) N — 50, and c) N = 55. The spots indicate the locations of the zeros of the canonical 
partition function. 



